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ABSTRACT 

The coalescence of compact binaries containing neutron stars or black holes is one of the most 
promising signals for advanced ground-based laser interferometer gravitational-wave detectors, with 
the first direct detections expected over the next few years. The rate of binary coalescences and the 
distribution of component masses is highly uncertain, and population synthesis models predict a wide 
range of plausible values. Poorly constrained parameters in population synthesis models correspond 
to poorly understood astrophysics at various stages in the evolution of massive binary stars, the pro¬ 
genitors of binary neutron star and binary black hole systems. These include effects such as supernova 
kick velocities, parameters governing the energetics of common envelope evolution and the strength 
of stellar winds. Observing multiple binary black hole systems through gravitational waves will allow 
us to infer details of the astrophysical mechanisms that lead to their formation. Here we simulate 
gravitational-wave observations from a series of population synthesis models including the effects of 
known selection biases, measurement errors and cosmology. We compare the predictions arising from 
different models and show that we will be able to distinguish between them with observations (or the 
lack of them) from the early runs of the advanced LIGO and Virgo detectors. This will allow us to 
narrow down the large parameter space for binary evolution models. 

Subject headings: binaries: close — gravitational waves — methods: data analysis — stars: black 
holes 


1. INTRODUCTION 


ometers are currently being commissioned an d should 
begin observing runs in 2015 ( Aasi et al.|2013b ) with the 
sensitivity increasing gradually over a number of years 
before reaching their design sensitivity near the end of 
the decade. These gravitational-wave (GW) observato¬ 
ries will be an order of magnitude more sensitive than the 
first generation observato ries and are expecte d to yield 
the first GW detections (Abadie et al. 2010) and her¬ 
ald the beginning of GW astronomy. In GW astronomy 
we are interested in the emission of gravitational radi¬ 
ation from astrophysical sources. One of the primary 
sources of GWs for aLIGO is the coalescence of compact 
binaries - binary neutron star (BNS), neutron star-black 
hole (NSBH) and binary black hole (BBH) systems. 

The orbits of these systems decay due to rad iation re¬ 
action (Peters & Mathews 1963 |Peters 1964), causing 
the two objects to spiral in towards one another. During 
the final orbits and merger, these sources emit a large 
amount of gravitational radiation, and this will be ob¬ 
servable by aLIGO and AdV. The gravitational wave¬ 
form emitted by the binary can be modelled with great 
accur acy using the post-Newtonian formalism (Blanche! 
2006). Gloser to merger, full numerical simulations are 
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required to tr ack the binary e v olution and ca l culate the 


The Advanced LIGO (aLIGO) (lAasi et al. 

2015 

) and 

Advanced Virgo (AdV) (|Acernese et al. |2U1 

5) second 

generation, kilometre-scale ground based laser interter- 


waveform (see Hannam (2009); Hinder (2010); Sperhake 


et al. (|2013 ) for overviews). By combining the insights 


of post-iNewtonian theory and numerical modelling, a 
large range of analytic/semi-analytic approximate wave¬ 
form mod els have been developed over the past few y ears 
(see, e.g., Buonanno et al.| (2009) and Ohme (2012) for 
an overview). These models now provide accurate wave- 
forms over a large fraction of the parameter space of non- 
precessing BBHs. In particular, they provide accurate 
waveforms for signals with a range of mass ratios and 
also c over the space of aligned spins. There is ongoing 


work (Hannam et al. 2014 Pan et al. 2014) to extend 


these to the full parameter space that incorporates spin- 
induced precession of the binary orbit. 

The availability of accurate waveform models makes 
a matched f ilter search of thes e signals feasible (Babak 


et al. 2013 

tract the piiysical parain eters of the binary system from 


Aasi et al. 2013c) and allows us to to ex- 


the observed GW signal (Aasi et aL]|2013a Veitch et al. 


2014). The observed sky location and orientation of the 


netic counterparts of GW systems I 

Abadie et al.||2012a 

Singer et al. 2014 

Aasi et al. 

2014 

Clark et al.l|2014). 


binary components will shed light upon the formation 
and evolution of the binary by comparing the observa¬ 
tions with predictions from stellar evolution models. We 
expect the majority of systems to be observed with rel¬ 
atively low signal-to-noise ratio (SNR) and consequently 
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the parameters will not be measured wit h great accuracy 
(Hannam et al.||2013 Ohme et al.||2013 ). For an individ- 
ual binary, the chirp mass ot the system — a combination 
of the two masses that determines the rate at which the 
binary evolves — can be measured with good accuracy 
(Cutler & Flanagan 1994 Hannam et al. 20131, while the 
mass ratio and spins are unlikely to be well constrained. 

In addition, there is significant uncertainty in the as- 
trophysical mass and spin distributions of black hole bi¬ 
naries. Thus, it seems unlikely that the measurement of 
parameters from individual systems will significantly im¬ 
pact our understanding of black hole binary formation. 
Instead, it will require the measurement of parameters 
from a population of signals to significantly constrain 
compact binary formation and evolution models. In this 
paper, we consider how this might be done and what 
we expect to learn with the observations from the early 
aLIGO and AdV runs. 

Compact binaries can be formed as a result of the 
evolution of isolated massive binaries (where the com¬ 
ponents have initial masses > 8Mq) or can be formed 
dynamically (i.e., in dense globular and nuclear star 
clusters) from binary-single star interactio ns between 


compact remnants an d primordial binaries (Mandel & 
O’Shaughnessy 20101. While the key stages of the bi¬ 


nary evolution are well understood, there are significant 
uncertainties in the details of the process. Population 
synthesis codes attempt to model these uncertainties us¬ 
ing empirical prescriptions. These models contain nu¬ 
merous parameters which are not well constrained relat¬ 
ing to astrophysics such as stellar winds, supernova kicks 
imparted on black holes at birth and common envelope 
binding energy among others. Varying these parameters 
will have a significant impact on both the predicted rate 
of compact binary mergers, as well as the distribution of 
expected masses and spins of the compact remnants that 
comprise the binary (Dominik et al.||2012). 

In this paper, we introduce a straightforward model se¬ 
lection method to distinguish between various formation 
and evolution scenarios. We focus on the two parame¬ 
ters that will be best measured: the overall rate of binary 
mergers and the chirp masses of the observed binaries. 

Furthermore, we restrict attention to BBHs as, based 
upon the recent population synthes is models, these are 
predicted to be the m ost numerous (Voss fc Tauris|200j 


Dominik et al.|[2012). We caution, however, that detec¬ 


tion rates are highly uncertain and previous papers have 


argued that there will be essentially no BBHs (Belczynski 
et al. 2007 Mennekens & Vanbeveren 20141. This triv¬ 


ially means that any detections ot merging BBHs will 
allow models predicting a dearth of such systems to be 
ruled out, shedding light on the astrophysical assump¬ 
tions made therein. Beyond that, we show how, in addi¬ 
tion to the merger rates, the broad range of BBH chirp 
masses predicted by population synthesis models encodes 
information about the BBH formation mechanisms. 

There have been many studies performed over the last 
decade that have made use of either one or both of these 
pieces of information t o distinguish between comp eting 


astrophysical models. Bulik & Belczyski (2003) used 


a Kolmogorov-Smirnov test to compare simulated GW 
chirp mass measurements to a series of predicted ob¬ 
served distributions from population synthesis models. 
They find they can distinguish many models with ^ 100 


observations, a find ing we confirm in the present study. 
Kelley et al. (2010| use a Bayesian approach introduced 


Mandel 12010[) to show how one can use GW observa- 


tions along with dark matter simulations to distinguish 
between different natal kick-velocity models, and again 
find they require 0(100) observations to distinguish be¬ 
tween models. 


Belczynski et al. (2012|) discuss using upper limits on 
binary merger rates to disti nguish between popu lation 


synthesis models. Recently, Mandel et al. 


(2015) have 


shown how one can use population synthesis models 
along with GW observations of binary mergers to mea¬ 
sure the relative rate of BNS, NSBH an d BBH merg¬ 
ers with 0(1 0) observations. In addition. Messenger & 
Veitch (2013) show how one should use all of the intorma- 
tion available to avoid selection biases when attempting 
to make inferences about distributions of rates and pa¬ 
rameters of merging binaries. 

More sophisticat ed techniques have also been discussed 


in the literature. O’Shaughnessy (2013) introduces a 


framework to incorporate measurements of both the 
merger rate and parameter distributions of GW obser¬ 
vations, and compares these to a set of population mod¬ 
els which sparsely sample the relevant parameter space. 
A sim ilar techni que is used in [Mandel fc 0’Shaughnessy| 
(|2010|) (see also|Mandel et al.| p010p). 

Here, we introduce a fast, simple method to make infer¬ 
ences about astrophysical models using information from 
GW observations. The method is general, and could be 
applied to any set of binary evolution models. We illus¬ 
trate its utility by evaluating the ability to di stinguish 


between a s uite of population synthesis models (Dominik 


et al. 2012). For concreteness, we restrict attention to 
the expected results from th e early observing r uns of the 


advanced GW detector era ( Aasi et ^|2013b ). 

Population synthesis models typically predict the 
galactic rate of binary mergers and the parameter distri¬ 
butions. From this, we model the observed distribution 
by accounting for observational bias: GW detectors are 
able to observe signals from higher mass systems to a 
greater distance. Additionally, we incorporate cosmolog¬ 
ical effects that lead to a red-shifting of both the observed 
masses and the observed merger rate. Finally, we model 
measurement errors and uncertainties inherent in the ex¬ 
traction of the signal from a noisy data stream. For each 
population synthesis model, we generate an expected ob¬ 
served rate and associated mass distribution. 

Based on simulated observational results, we can use 
model selection to differentiate between the various mod¬ 
els. To give a sense of what we can expect, we simu¬ 
late results from the early aLIGO and AdV observational 
runs. To do this, we choose one of models from a suite of 
population synthesis models to play the role of the uni¬ 
verse, and draw GW observations of BBHs from it, ac¬ 
counting for known observational biases and anticipated 
measurement errors. We then compare these observa¬ 
tions to the full suite of population synthesis models and, 
starting with a uniform prior on the models, we compute 
the posterior probability for each model. 

While the results that we present are limited to these 
specific scenarios, the method we introduce is general 
and could easily be applied to the predictions from any 
population synthesis model and the results (predicted or 
observed) from any detector network. We also caution 
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the reader that the models of|Dominik et al. (2012 1 rep¬ 
resent the most optimistic predictions ot h5BH nierger 
rates, with all models predicting a detection within the 
first two aLIGO and AdV science runs. Lower merger 
rates would lead to observations of BBH mergers only 
in later runs at, or close to the design sensitivity of the 
detectors. Fo r an overview of rate predictions for aLIGO 
and AdV see |Abadie et al.| (|2010|). 

This paper is structured in the following way. In Sec- 
tion[^we give a brief review of compact binar y for mation, 
and introduce the models we use in Section l2.2l In Sec- 
tionj^we describe our algorithm for accounting for known 
selection biases, converting an intrinsic chirp mass dis¬ 
tribution to a predicted observed distribution. Section 
shows how to use information from the two well mea¬ 
sured parameters — the chirp mass and the merger rate 
— to distinguish between population synthesis models. 
In Sections [^&[^ we show what we may be able to learn 
about binary evolution using GW observations of binary 
black holes from the first two aLIGO and AdV science 
runs. Finally in Section[^we conclude and suggest areas 
which require further investigation. 

2 . GOMPAGT BINARY FORMATION AND 
EVOLUTION 

In this section, we provide a brief review of isolated bi¬ 
nary evolution, highlighting the poorly understood stages 
of the evolution, which lead to the uncertainties in the 
predicted merger rates and mass distributions of t he bi¬ 
naries. For more i nform ation see a review such as IPost-l 
nov & Yungelsonl (|2006|). 


2.1. General overview 

For a single star, its evolution is solely determined by 
the zero-age main sequence (ZAMS) mass and composi¬ 
tion. However, the majority of massive stars exist in bi¬ 
naries or multiple systems, with > 70% of massive 0-type 
stars exchanging mass with a companion during their life¬ 
time (Sana et al.| 2013 Duquennoy & Mayor 19911. In 
this case, the evolution is no longer straightforward, and 
can lead to a plethora of exotic systems. Here we give 
one possible evolutionary pathway for a massive binary; 
many alternati ve pathways also exis t (see for example 
Tables 4 & 5 in Dominik et al. (2012) for a summary). 

Gonsider a binary in which both stars have ZAMS 
masses > 8Mq . The initially more massive star (the pri¬ 
mary) in the binary will evolve off of the main sequence 
first since it has the shorter lifetime. As it evolves, its 
radius expands until it fills its Roche Lobe as a giant and 
begins to transfer mass to the companion (the secondary) 
star, stripping the primary’s hydrogen outer layers and 
leaving a He/Wolf-Rayet star. Already the evolution of 
the binary is different to that of single stars since the 
companion can change its mass considerably, leading in 
some cases to a reversal of the mass ratio. If the core 
is massive enough, the primary will then collapse in a 
supernova, and leave behind a compact remnant — ei¬ 
ther a neutron star or a black hole depending on the 
pre-supernova core mass. 

In stellar evolution models, the distinction between col¬ 
lapse to a neutron star or a black hole is made via mass 
alone, with the maximum allowed mass of a neutron star 
being one of the free parameters. In reality, the maxi¬ 
mum neutron star mass is set by the unknown neutron 


star equation-of-state. The maximum observed neutron 


stars have masses aro und 2 M 0 (Demorest et al. 20I0_ 
Antoniadis et al.|2013). Gausality and General relativity 


require the maximum neu tron star mass to be < 3 . 2 M 0 
( Rhoades &: Ruffmi||1974 ). 

The mechanism of the supernova itself is intensely 
studied but still not fully understood. If the supernova 
is asymmetric (due to asymmetric mass loss or neutrino 
emission) the resulting neutron star can be given a na¬ 
tal kick velocity due to the conservation of momentum, 
which is of the order 25 0 km s~^ for galactic neutron 
stars (Hobbs et al. 2005). It is unclear whether black 
holes also receive a kick of this magnitude or whether 
mass falling back onto the black hole reduces the size o f 
this kick sign ificantly (see for e.g. Repetto et al. (2012); 


Janka (2013|)). 

It the system survives the first kick, then the secondary 
will begin to evolve. The compact remnant accretes mat¬ 
ter from the stellar wind of its companion, becoming a 
luminous X-ray source. At this stage, the binary may 
be observable electromagnetically as a high-mass X-ray 
binar y. Although the th eory of stellar winds is fairly ro¬ 
bust ( Castor et al.||1975 ), the str ength of stellar winds i n 
these systems remains uncertain ( Lepine fc Moffat|2008 ). 

As the secondary continues to evolve, it will continue 
to expand and fill its Roche Lobe. If the mass transfer 
through Roc he Lobe Overflow is unstable, a comm on en¬ 
velope phase (Ivanova et al. 2013 Paczynski 1976) can be 
initiated. This is where both the compact remnant and 
the core of the secondary orbit within the secondary’s hy¬ 
drogen outer layers. The common envelope is the least 
well understood phase in the evolution of binaries. The 
common envelope is usually parametrised in one of two 
fashions; the a prescription ( Webbini^|1984 1 f ocusing on 


conservatio n of energy, or the 7 prescription (Nelemans 
et al.|[ 2000 ) focusing on conservation of angular momen¬ 


tum. i'he core and compact object spiral in towards one 
another on a dynamical timescale due to drag, and or¬ 
bital energy is used to eject the envelope. This stage is 
responsible for dramatically reducing the orbital separa¬ 
tion in the binary. 

If the binary survives the common envelope, the core of 
the secondary can then go supernova, potentially impart¬ 
ing a second kick on the system (although it is generally 
less likely to unbind the system since the orbital veloc¬ 
ities are now much higher). Finally, a compact binary 
remains containing neutrons stars and/or black holes. It 
is these systems which then inspiral towards one another 
and merge due to radiation reaction, and will be observed 
in GWs by aLIGO and AdV. 


2.2. Detailed binary evolution models 

Population synthesis codes are Monte-Garlo simula¬ 
tions that evolve large ensembles of primordial binaries 
via semi-analytical prescriptions, taking as input param¬ 
eters corresponding to the poorly understood astrophys- 
ical stages outlined above. Binary population synthesis 
models can be used to try to understand the effects of 
these uncertainties on binary evolution, and on the re¬ 
sultant population of compact binaries. One way to ex¬ 
ploit the information contained in GW observations of 
coalescing BBHs is therefore to compare the measured 
properties of a population to population synthesis model 
predictions. 
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Table 1 

Summary of population synthesis models. 


Model 


Physical difference 


Standard 


Variation 1 
Variation 2 
Variation 3 
Variation 4 
Variation 5 
Variation 6 
Variation 7 
Variation 8 
Variation 9 
Variation 10 
Variation 11 


Maximum neutro n star mass = 2.5 Mq, rapid 
supernova engine ||Fryer et al.|2012|| , physically 
motivated envelope binding energy IXu fc Li| 
20101, standard kicks a = 265 km s^^ 


Very high, fixed envelope binding energy^ 
High, fixed envelope binding energy^ 

Low, fixed envelope binding energy^ 

Very low, fixed envelope binding energy^ 
Maximum neutron star mass = 3.0 Mq 
M aximum neutron mass = 2.0 Mq 
Reduced kicks a = 123.5 km s“^ 

High black hole kicks, /^ = 0 
No black hole kicks, /{, = 1 
Delayed supernova engine (|Fryer et al.|2012|l 
Reduced stellar winds by factor of z 


Note. — Models presented in|Dominik et al.| l|2012||, with 
parameter variations indicated in ttie second column which 
broadly relate to the uncertainties in binary evolution discussed 
in the text. All other parameters retain their standard model 

value. _ 

^ See Section 2.3 in |Dominik et ar] ( |2012[ | for details 

For this study we use a set of publicly availablFI pop¬ 
ulation synthesis models presented in |Dominik et al. 
(2012), pro duced using the StarT rack population syn¬ 
thesis code (Belczynski et al. 2008). Predicted chirp mass 
distributions and merger rates ot BNS, NSBH and BBH 
systems are provided for a range of input physics. 

The relative rates of BNS, NSBH and BBH merg¬ 
ers are uncertain. Although BBH systems are more 
massive (and consequently detectable to a greater dis¬ 
tance), much more massive stars are needed in order to 
form them, and the initial-mass-function (IMF) falls off 
rapidly at high masses, meaning these stars are rarer. 
It is also worth noting that no BBH has ever been ob¬ 
served, although system s which may be progen itors for 
them such as Cyg X-3 (Belczynski et al. 2013 ), IC 10 
X-1 (Bulik et al. 2011) and JNGC 300 A-1 (Crowther 
. |2010|) have been studied and provide some limits 
on BBH merger rates. The population synthesis model 
we are utilising predicts that BBH detection rates will 
dominate over those for BNS and NSBH. Based on this, 
and to keep the analysis simple, we restrict our attention 
to BBH systems. It would be relatively straightforward 
to extend the framework we introduce to include all GW 
observations of binary mergers. 

We use the set of 12 population synthesis models for 
which predicted rates and mass distributions are avail¬ 
able. These models are summarised in Table H The 
standard model assumes a maximum neutron star mass 
of 2.5 MfT>, uses th e rapid supernova engine detailed in 


Fryer et al. (2012), physically motiv ated common enve- 
lope binding energy (Xu & Li 2010), and kick velocities 
for supernova remnants drawii frorri a Maxwell distribu¬ 
tion with a characteristic velocity of cr = 265 km s“^. 
There are then eleven variations, in each of which one 
of the above parameters is varied: the first four varia¬ 
tions consider changes in the energetics of the common 
envelope phase, the next two vary the maximum mass of 
neutron stars, three more change the kick imparted on 

^ http: //www. syntheticuniverse. org 


the components during collapse to a neutron star or black 
hole and the final two consider a delayed supernova en¬ 
gine and a change in the strength of stellar wi nds. The 
models are described in detail in section 2 of IDominikI 


et al. (2012). 

We expect that in general, the true universe will not 
match one of a small set of models, but will lie in between 
these models (or potentially outside of them if additional 
unmodelled physics is required to accurately describe bi¬ 
nary evolution). Assumptions that are not varied in these 
models, but which may have a large impact on the re¬ 
sultant BBH distribution include distributions of the pa- 
rameters of primordial binar ies (IMF, orbital elements 
de Mink fc Belczynski||2015 )), tides, stellar rotation (de 
dink et al. 2(J13() and magnetic fields. Here we neglect 
these additional considerations and investigate how one 
can differentiate between a small suite of population syn¬ 
thesis models using GW observations of BBHs. A full 
treatment of these additional properties has the potential 
to significantly impact stellar evolution models and may 
well lead to degeneracies whereby significantly different 
astrophysical models predict comparable populations of 
binaries. 

Since calculating population synthesis models can be 
computationally expensive, the models are discretely 
sampled over a large range of parameter space (in some 
cases orders of magnitude) in an attempt to bracket the 
truth. Furthermore, each of the models used in this study 
varies only one parameter from its standard value at a 
time. It is quite likely that the true values of many of 
these pararn e ters w ill differ from those presented in |Do-| 
minik et al. (2012), resulting in a population that does 
not match any ot the ones included here. Varying combi¬ 
nations of parameters will also need to be studied, as this 
may lead to issues with degeneracies in which combina¬ 
tions of parameters can be determined from GW observa¬ 
tions. To be able to reliably extract the details of stellar 
evolution from GW observations, one would require to 
have models calculated on a dense enou gh grid that one 


can perform interpolation betwee n them (O’Shaughnessy 
2013 lO’Shaughnessy et al.||2008 ). 


2.2.1. Metallicity 

Each model is calculated at solar {Z — 0.02 = ^ 0 ) and 
sub-solar {Z — 0.002 = O.IZq) metallicities. In addition, 
there are two submodels that differ in the way binaries 
entering into a common envelope when one of the stars is 
on the Hertzsprung gap are handled (see section 
We choose to use a 50-50 mixture of the 


so 


2 . 2 . 2 ). 


ar and 


sub-solar models as used in Belczynski et al. (2010), mo¬ 
tivated by results from th e Sloan Uigitai Sky Survey 
(SDSS, ( jPanter et al.|2008 )) showing that star formation 
is approximately bimodal with half of the stars forming 
with Z ~ ^0 and the other half forming with Z ~ O.IZq. 
For the future, it would be desirable to include a more 
thorough treatment of the metallicity distribution, in¬ 
cluding its evolution with cosmic star formation history 


as done in Dominik et al. (2013 2014). 


Although metallicity in the local universe may be bi¬ 
modal, one still expects a smooth distribution of metal¬ 
licities to exist. Using only a discrete mixture of so¬ 
lar and sub-solar metallicity predictions may give rise to 
non-physical peaks or sharp features in the chirp mass 
distributions which may artificially aid in distinguishing 
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between them (Dominik et al.||2014). However, in prac¬ 
tise we find that these peaks are suffi cient ly smoothed 
out by measurement errors (see section 3.4). 

Studies have shown that the majority oFHBHs observ¬ 
able by aLIGO were formed within ^ 1 Gyr of the Big 


Bang (jDominik et al. 2013 20141, when the metallicity of 
the uniyerse was distinctly low er. This is due to a num - 
ber of reasons (see for example Belczynski et al. (2010)). 
It is easier for supernova progenitor stars to remain mas¬ 
sive at lower metallicities due to weaker stellar winds 
compared to at solar metallicity. Also, many potential 
BBH progenitor systems merge prematurely at higher 
metallicities during the CE phase since the secondary is 
likely to be on the Hertzsprung Gap, whereas at lower 
metallicities the secondary does not expand enough to 
initi ate a GE even t until it is a core-helium burning star 
(2000| for the effect of metallicity on 
nese BBHs are formed with long delay 


(see Hurley et al. 
stellar radius), i'. 


times such that they are only merging now. One there¬ 
fore needs to include the time evolution of metallicity 
to accur ately model the expe cted population of BBHs 
mergers (Dominik et al.|[2013). 


2.2.2. Fate of Hertzsprung Gap donors 

The Hertzsprung gap is a short lived (Kelvin- 
Helmholtz timescale) stage of stellar evolution where a 
star evolves at approximately constant luminosity across 
the Hertzsprung-Russell diagram after core hydrogen 
burning has been depleted but before hydrogen shell 
burning commences. 

Whilst on the main sequence, stars are core burning 
hydrogen, and do not possess a core-envelope separation 
as the helium core is still being formed. Therefore, if 
a main sequence star enters into a common envelope, 
orbital energy is dissipated into the whole star, rather 
than just the envelope, making ejecting the envelope ex¬ 
tremely difficult. It is therefore expected that any star 
entering into a common envelope phase whilst on the 
main sequence will result in the two stars merging pre¬ 
maturely in an event which is not interesting from a GW 
astronomy point of view. 

Eor stars that are on the Hertzsprung gap, the situation 
is not so clear. The helium core begins contracting whilst 
the envelope of the star expands. It is unclear if there 
is sufficient core-envelope separation on the Hertzsprung 
gap for a star entering a common envelope phase to have 
its envelope ejected, or whether it would suffer a similar 
fate to a main sequence star. 

The fate of Hertzsprung Gap donor s is another of 
the un certainties that is investigated by |Dominik et al. 
(2012). In t he optimistic submode l (referred to as sub- 


mode 


A ir 


Dominik et al. (2012[)), the authors ignore 


the issue a nd calculate the common envelope energetics 


as normal (|Webbink|| 1984|). In the pessimistic submodel 
(referred to as subrnodel B), any binary in which the 
donor is on the Hertzspung gap when the binary enters 
into a common envelope phase is assumed to merge. This 
tends to reduce the number of merging binaries (and thus 
the rates) compared to the optimistic model. It is un¬ 
likely that either of these models is accurate, as the fate 
of a Hertzsprung gap donor will depend on the inter¬ 
nal structure of the star as it enters the common enve¬ 
lope phase. Nonetheless, submodels A and B provide 
upper and lower bounds, respectively, on the number of 


Hertzsprung gap donors forming BBH. 

In this paper, we compare results for the twelve models 
listed above for both the optimistic (submodel A) and 
pessimistic (submodel B) Hertzsprung gap evolution. 


3. PREDIGTED OBSERVED DISTRIBUTIONS 

Eor each of the models described above, we are given an 
expected rate of binary mergers per Milky Way Equiva¬ 
lent Galaxy (MWEG), as well as a distribution of binary 
parameters (notably the component masses). The popu¬ 
lation of BBHs observed by the advanced GW detectors 
will differ from this underlying intrinsic distribution due 
to the following observational effects. 

(a) The GW signal from binaries at large distances will 
be red-shifted due to the expansion of the universe 
which consequently leads to a shifted measurement 
of the binary’s total mass. 

(b) The GW amplitude scales with the binary’s total 
mass, thus binaries with heavier components will be 
observable to greater distances, provided their signal 
still lies in the sensitive frequency region of the detec¬ 
tor, which leads to an increased number of observed 
high-mass systems. 

(c) Due to the presence of noise in the detector the best- 
measured parameters will differ from the binary’s 
intrinsic parameters which effectively blurrs the ob¬ 
served distribution. 


We take all three effects into account and calculate the 
distribution of parameter we expect to observe. Our ap¬ 
proach is con sistent with pre v ious m ethods in the lit¬ 


erature (e.g., Dominik et al. (2014)), apart from how 


we account for measurement errors across the parameter 
space. Eor completeness, in the remainder of the section, 
we briefly recap how these effects are accounted for and 
the observed distribution obtained. 


3.1. Detectability 

We model the GW signals by the dominant harmonic 
only, which is sufficient for the majority of black hole 
systems we are considering (|Gapano et al.|2014| |Bustillo 
et al.||2015|). The sign al observed m a GW def ector can 
then be exjaressed as ( [Eairhurst k, Brad^|2008[ ) 

h{t) = ^ [/io(t)cos$-h/i^/2(t)sin$] , (1) 

n'eff 


where is called the effective distance, $ is the coa¬ 
lescence phase as observed in the detector and are 

the two phases of the waveform which are offset by 7r/2 
relative to each other [equivalently, their Eourier trans¬ 
forms obey ho(/) = *^ 7 r/ 2 (/)]- The effective distance is 
defined as 

DeS = , ^ ^=- ( 2 ) 

y E+(l + cos^ tY cos^ L 

Dl is the luminosity distance to the binary, are 

the detector response functions and t is the binary in¬ 
clination angle. The maximal (and circularly polarised) 
GW signal is obtained when the signal is directly over¬ 
head the detector F+ = 1; Ex =0 and with (, = 0 , tt 
corresponding to a face on signal. 
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The effective distance is in versely proportional to the 


SNR, p, which is de fined as ( [Cutler fc Flanagan||1994 
Poisson fc WilljfTg^ l: 


p^=A 




\HfW 

Snif) 


df, 


( 3 ) 


where h{f) is the frequency-domain gravitational wave¬ 
form and the detector noise power spectral density is 
denoted by S'„(/). We choose a lower cutoff frequency of 
/low = 20 Hz, suitable for the early advanced detectors. 
The SNR at which a signal can be detected will depend 
upon the details of the detector network, including the 
sensitivities of the detectors as well as the character of 
the data — non-stationarities in the data make it more 
difficult to distinguish candidate signals from the back¬ 
ground noise. However, for studies such as this, it is con¬ 
venient to choose an approximate threshold. Experience 
has shown that a network SNR of 12 is a pproximately 
where we might expect to make a detection (jAbadie et al. 


2012b Aasi et al. 2013b). This corresponds to an SINR 
of around S m each of tne LIGO detectors in the early 
science run^ For the studies presented in this paper, 
we use this simple, single detector threshold to decide 
whether a signal would be observed by the detector net¬ 
work. 

Given a model for the waveform, h{t), we can calculate 
the maximum effective distance to which the signal could 
be detected. This is known as the horizon distance, Dh, 
and corresponds to the maximal distance at which the 
signal could be observed if it is optimally oriented and 
overhead. To calculate the horizon distance we us e the 
phenomenolog i cal wa veform model introduced by jSan-j 


tamarfa et al. ( 2010 ) that includes the inspiral, merger 
and ringdown sections of the waveform calibrated using 
numerical relativity. The model provides the waveform 
h(f) in the frequency domain as a function of the bi¬ 
nary’s total mass M, its symmetric mass ratio rj and an 
effective total spin parameter, y. 

The mass parameters of the binaries are characterized 
in terms of the best measured parameter combination, 
the so called chirp mass M, which is a combination of 
the component masses mi and m 2 . 


M = MX’'*. 

(mi -I- m 2 ) 


( 4 ) 


where M = mi + m 2 is the total mass, and 77 is the 
symmetric mass ratio, 


T] = 


mim2 


[mi -\- m2)^ 


< 0.25. 


( 5 ) 


For an equal mass binary mi = m 2 = m, the chirp mass 
M « 0.87m. In the remainder of the paper, we will focus 
on the predicted and observed chirp mass distributions, 
and not consider mass ratio or spin. 

Our aim is to predict the observed chirp mass distri¬ 
bution given the intrinsic model prediction, and compare 
these with observations. 


^ For the early science runs, we expect the LIGO detectors to be 
about twice as sensitive as Virgo so, on average, one might expect 
a threshold event to have SNR of 8 in each of the LIGO detectors 
and 4 in Virgo 



chirp mass [Mg] 

Figure 1. Horizon distance in Gpc for nonspinning BBHs as a 
function of chirp mass and symmetric mass ratio assuming a single 
detector with the early aLIGO noise spectrum. 


Throughout most of this paper, we assume the 


aLIGO (circa 2015) noise spectrum (Aasi et al._ 
2013b|) representing the expected sensitivity ol a 


early 


2010 

:iuD 

during its first observing runs. A plot of the horizon dis¬ 
tance as a function of the chirp mass and the symmetric 
mass ratio is given in Figure It encodes the farthest 
distance to which a BBH with the given parameters can 
be seen. The horizon dista nce will also be a fun ction of 
the black hole spins. Since Dominik et al. (2012) do not 
provide individual spin information in their catalogues, 
we set the spin parameter to zero for simplicity when sim¬ 
ulating signals in our synthetic universe. Our ignorance 
of the spins may lead to systematic biases, as high spins 


can n oticeably affect the horizon distance (Ajith et al. 


2011 ) and could chang e the rate of observed signals by a 


factor of two or three (Dominik et al. 2014). One could 
incorporate this lack of knowledge by assuming a spin 
distribution for black holes and margnializing the result 
over the spins. We defer this to a later study when more 
informed spin priors (observationally motivated or from 
population synthesis) can be incorporated. 

Not every binary within the horizon will be detected, 
as DeB is location and orientation dependent. Under the 
assumption of a uniform distribution over the sky and a 
uniform source orientation, however, we can numerically 
calculate the fraction P(^) of systems with 


+ cos^ i)^/4 -I- Fx cos^ 6 > (6) 

with ^ G [0,1]. Note that P(^), which we can inter¬ 
pret as a cumulative distribution function, is independent 
of the binary’s masses, and we will use it to determine 
what fraction of signals at a given luminosity distance 
is detectable, i.e., has an SNR larger than the detection 
threshold. 


3.2. Cosmological effects 

We simulate an expanding universe with sources dis¬ 
tributed uniformly and isotropically in comoving volume, 
which on scales of hundreds of Mpc is a valid assumption. 
Since the frequencies of any signal become increasingly 
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redshifted with growing distance between source and de¬ 
tector, the total chirp mass measured at the detector is 
shifted by 

M* =Mil + z), (7) 


where z denotes the redshift. Assuming zero curvature 


and s tandard cosmological parameters (Bennett et al. 


20141 


flM = 0.286, OA = l-nM, Ho = 69.6, (8) 

we calcu late the com oving distance as a function of the 
redshift ( Hogg|[T9M ), 


Dc{z) = 


dz' 


H 


0 Jo 


•\/f^M(l 


,/'|3 


r^A 


( 9 ) 


Here, c denotes the speed of light. 

The catalogues bylDominik et al. (2012) provide large 
sets of binaries characterized by their intrinsic chirp mass 
A4 and symmetric mass ratio ly. When we distribute 
them uniformly in comoving volume, the observed chirp 
masses M* are redshifted according to ([^. This implies 
that the maximal distance to which they can be detected 
changes as it is the observed parameters, not the intrin¬ 
sic parameters, that determines the appropriate horizon 
distance. Since 


Dl = Dc{1 + z), (10) 

the maximal observable comoving distance satisfies 

D^--{M,7i,z){1 + z)=Dh{M*,v), ( 11 ) 

which we solve numerically for z. Note that the leading- 
order inspiral horizon distance behaves as Dh{M*) ~ 
hence 

7Vf5/6 

^ < ^h{M). (12) 

While the derivation of ( [T^ is only valid for low-mass 
systems, we find that is generally less than the 

static. Euclidean universe equivalent Dh- 


3.3. Detection rate and distance distribution 

We now assume binaries of a fixed model, distributed 
isotropically and uniformly in comoving volume, that 
merge at a constant comoving merger rate dens ity TZ (in 
MWEG~^ M yr~^) as given in the data sets by 


Dominik 


et al. (2012). To convert these numbers into a detection 
rate tor aLlGO, we proceed as follows: 

First, the comoving merger rate TZ per MWEG has to 
be multiplied by an average gal axy density which we tak e 
as Pa ~ 0.0116 Mpc“^ following Kopparapu et al. (2008). 


Next, we must calculate the effective volume in which 
each binary is observable, by integrating the number of 
observable binaries as a function of Dc- As the distance 
increases, the area of the corresponding sphere increases 
as but the fraction of binaries that are oriented such 
that their signal is sufficiently loud for detection (that 
is, Heff < Dh) becomes smaller. Finally, due to the 
redshifted time, 

tL =tcil + z) (13) 


the observed merger rate for binaries at redshift z > 0 is 
less than the comoving merger rate. Thus, the effective 


redshift 



Figure 2. The probability distribution in comoving distance for 
detectable BBHs with mi = m 2 = IOMq. The solid (blue online) 
curve takes cosmological effects into account (see text) while the 
dashed line assumes a static, Euclidean universe (i.e. local universe 
approximation). 


volume for a binary with chirp mass Ai is 


Ves{M) = 47r 






l-fz \DH{M*,r]) 


dDc, 


(14) 

where D^^^ is defined by (|11|). The function P, intro¬ 
duced in Eq. (1^, gives the Traction of suitably oriented 
binaries (i.e., those giving an SNR greater than 8) and 
(1 -I- z)~^ accounts for the difference between apparent 
and comovin g m erger rate density. We note that the 
integrand in (14) can be interpreted (up to a normali¬ 
sation) as the observed distance distribution for binaries 
with fixed intrinsic parameters. 

The average detection rate for each model is given by 


N = TZ X pc y. 14ff, 


(15) 


where VeS denotes the average effective volume, with the 
average taken ov er all binaries in a given mo del. We take 
TZ and pn f rom Dominik et al. (2012) and Kopparapu 
et al.| (|2008|), respectively. 

Figurej^shows this distribution for an equal-mass BBH 
with TOi = 7712 = IOMq. For comparison, we include the 
equivalent curve for a static, Euclidean universe, where 
Dl = Dc and (14) simplifies to the case z = 0. As 


expected, both curves agree for low redshift, but as we 
have noted above, there are fewer binaries seen at large 
comoving distances if the expansion of the universe is 
taken into account. This effect becomes increasingly im¬ 
portant for larger distances, i.e., for high-mass binaries 
and more sensitive detector configurations. 

The effective volume in which binaries with fixed pa¬ 
rameters are detectable changes considerably across the 
parameter space. This leads to an observational bias in 
favour of systems with large volume reach. We incor¬ 
porate this effect by re-weighting the chirp mass distri¬ 
bution of binaries acc ording to their i ndivid ual effective 
volumes. In practice, Dominik et al.j(2012) provide the 
data for each of their models in form of a discrete set of 
binary parameters. For each of those binaries, we cal¬ 
culate the integer part of I4ff/I4g“ and add that many 
copies of the binary to our new set of observable param¬ 
eters. Here, denotes the smallest effective volume 
across all binaries in the set, and only one copy of the 







































binary with this smallest effective volume is keptj^ 

Finally, for each binary in our new set, we draw a co¬ 
moving distance from the distribution underlying 14ff. 
From this distance, we then infer the redshift and change 
from A4 to the observable redshifted chirp mass A4* ac¬ 
cording to 0. Our discrete representation of observable 
binaries then consists of multiple copies of the same in¬ 
trinsic systems, each however with a unique redshifted 
chirp mass. 

Note that an equivalent, but computationally more ex¬ 
pensive, procedure would be to randomly draw binaries 
from the intrinsic distribution, then draw comoving dis¬ 
tances and orientations for each binary within the total 
sensitive volume for the respective model and only keep 
those binaries with a detectable GW signal. Our method 
instead avoids disregarding any randomly drawn sources 
by drawing from the appropriate (distance/orientation) 
distribution of detectable signals. 


3.4. Estimating and Including Measurement Errors 


Including the observational bias discussed in Sec. |3.3 
in the chirp mass distribution still does not yield the 
distribution that one would expect to observe, because 
there will be a measurement error associated with each 
of the observations. Previous publications have mainly 
discussed a full Bayesian framework to combine mul- 


taint ies 

Mandel 2010 

Mandel & O’Shaughnessy 

2010l 

U’Shaughnessy 

2m 

). We instead assume a statistical 


value as detailed below. 

The accuracy of the parameters recovered during GW 
searches is limited by two factors. First, since we match 
to templates of the signals, the accuracy of recovered pa¬ 
rameters will be limited by the accuracy of the waveform 
models that used in the search. Second, the accuracy 
will be affected by statistical fluctuations of the noise 
in the measurement process. While we leave the fo rmer 
for dedicated stu dies such as|Buonanno et al.|(|2009|) and 
Nitz et al. (20131, we can estimate the uncertainty due to 


the latter using the well-known Fisher matrix estimate. 

Fisher matrix analyses rely on a linear approximation 
of signal variations and are valid for large SNRs. Nei¬ 
ther of the two assumptions is typically valid in realistic 


cations of violating these assumptions (Vallisnerif 

2008 

Rodriguez et al. 

2013 

Mandel et al. 

2014|). Here, 

how- 


method to distinguish BBH populations with GW ob¬ 
servations, we take Fisher-matrix predictions as a proxy 
for the width of posterior distributions of parameters ob¬ 
tained via a fully Bayesian analysis of the kind that wil l 
be performed on actual GW events (Veitch et al.||2014). 


In performing a population study ot the kind we per¬ 
form here, one should include not only a point estimate 
for parameters such as the chirp mass, but the full pos¬ 
terior from these parameter estimation routines. These 
posteriors c an then be com bined in the correct way, as 


described in Mandel (2010). The method we use here is 


^ We could keep more copies of the binary with the smallest 
effective volume and multiply the number of every other binary in 
the set accordingly, but tests showed that this has no effect on our 
results. 


essentially the point estimate approximation to the full 
analysis. 

We employ the same i nspiral-merger-ringdown model 
(Santamaria et al. 2010) for our Fisher-matrix calcula- 
tions as we used to simulate GW signals. We only con¬ 
sider variations of the intrinsic parameters: masses, time, 
phase and a model-specific single effective spin. We as¬ 
sume that these are also the parameters that are re¬ 
covered, at least initially by the GW search algorithm 
(see, e.g., the recently proposed se arch algorithm for 
nonpr ecessing, spinninng binaries by |Dal Canton et ah] 
(2014)). This assumption is likely to make our error es- 


timates too large since actual GW events will be fol- 
lowe d up by complex par ameter estimation routines (see 


space ot processing 

binaries (iVitale et al. |2014 

Ghatzi- 

ioannou et al. 

2015 

U'Shaughnessy et al.| 


). How- 


ever, since we only need an approximate error estimate 
that can be obtained in a fast and reliable way across 
the BBH parameter space, we choose to use the Fisher- 
matrix method here for nonprecessing binaries, and we 
neglect small correlations with extrinsic parameters such 
as sky location, orientation or distance. 

The characteristic standard d eviations in the measu re- 
ment process are estimated by (Poisson fc Will|[T995|) 


A0* = V(r-i)*,, 

dh 


(16) 

(17) 


where F^ is the Fisher information matrix and h = h (O^) 
is the waveform model. The inner product used in (17) 
is given by 


{g\h) = 4Re 


'/lo 


Snif) 


df 


(18) 


which is consistent with the SNR definition in (1^. The 
form of the waveform model we use allows us to cmculate 
the partial derivatives used in the definition of F^ ana¬ 
lytically, and we ensure numerical errors in the matrix 
inversion remain small [3 

The only parameter we use to distinguish BBH pop¬ 
ulations in this study is the observable, redshifted chirp 
mass, M*. The data sets of expected observable chirp 
masses that we p repared following the algorithm intro¬ 
duced in Sec. |3.3| shall now be skewed further by adding 
measurement errors to each binary in the data set. We do 
so by assuming a Gaussian distribution centred around 
the chirp mass value of each binary with a sta nda rd de¬ 
viation given by the Fisher matrix estimate (16). We 
evaluate the Fisher matrix at the appropriate observed 
chirp mass and mass ratio of the binary, setting the value 
of the black hole spins to zero (although we allow the 
spins to vary when calculating the Fisher matrix). This 
has a negligible effect on our results as the measurement 
accuracy for the chirp mass is only weakly dependent 
on the spins (Ohme et al. 2013). We randomly draw a 
sample from this distribution to re-define the measured 
chirp mass. Similarly, when we later simulate the uni- 


^ In fact, we find that no element of and F~^F deviates 

from the respective element of the identity matrix by more than 
10“^, in most cases the deviation is much less than this. 
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10 20 30 40 50 60 

chirp mass [Mg] 


Figure 3. Expected relative measurement errors in the chirp mass 
for an early configuration of aLIGO, SNR 10, calcul ated using the 
PhenomC waveform model ([Santamarfa et al.|2010|. 


verse with a particular model, each observation is drawn 
from the distribution that incorporates observational bi¬ 
ases, but the actually measured chirp mass is additionally 
offset following the Gaussian distribution that simulates 
measurement errors. 

The Fisher-matrix estimates scale inversely with the 
SNR, so we only calculate them once across the parame¬ 
ter space and scale them for each binary in the data set 
according to its SNR, which in turn is inferred from the 
distance and a randomly chosen orientation. Figure 
shows the chirp mass uncertainty at a constant SNR m 
10 across the parameter space for the early configuration 
of aLIGO. 

Figure illustrates the tra nsition from th e intri nsic 
BBH popmation, predicted by Dominik et al. (2012) for 
each of their models, to the expected observed chirp 
mass distribution. The main effect of the observational 
bias detailed in Sec. 13.31 is that the distribution becomes 
skewed towards high-mass binaries, and its support ex¬ 
tends to larger chirp masses due to the redshift of distant 
sources. The addition of measurement errors hardly af¬ 
fects the distribution at low chirp masses, simply because 
the errors are small compared to the typical variation of 
the distribution in this regime. For heavy systems, on the 
other hand, noise fluctuations introduce a non-vanishing 
chance of measuring chirp mass values greater than the 
largest (redshifted) chirp mass in each model. Hence, the 
main effect of introducing measurement errors is that the 
expected observed distributions show a characteristic tail 
at high chirp masses. 


4. GOMBINING MEASURED RATES AND GHIRP 
MASSES 

Given a set of BBH observations, for each model varia¬ 
tion U, we wish to calculate the posterior probability for 
that being the correct model. The information we gather 
about the correct model is twofold. First, we obtain a set 
of observed chirp masses {Al}, and second, we measure 
the rate of BBH detections by observing n binaries in 
a given observation period. In reality each observation 
will include measurements of additional physical param¬ 
eters of the system, such as the component spins (see 


Gerosa et al. (2013 2014) for information on how mea- 
surements of spin niisalignments can help to constrain 
astrophysical formation scenarios.) Including additional 
information from these other dimensions should help in 
distinguishing astrophysical formation scenarios. 

We simulate the observed population by choosing one 
of the model variations, adjusted to account for selection 
effects as described above, to describe the universe. We 
then draw n individual chirp mass measurements from 
this model, which form the data {Al}. The number of 
observations we assume is itself drawn from a Poisson 
distribution with a mean value that is dictated by the 
observation time and the merger rate of the model vari¬ 
ation we have selected to simulate the universe. 

With these measurements, {Al} and n, the posterior 
probability that model U is the correct model reads 


P(U|{Al},n) = 


P({Al},n|U)P(U) 
P({Al},n) ^ 


(19) 


where we have used Bayes’ Theorem. P(Ui) is the prior 
probability on model V), U({Al},n|Ui) is the likelihood 
of making these chirp mass measurements and measur¬ 
ing this detection rate given model variation Vj, and 
P({Al},n) is a normalisation factor called the evidence. 

Assuming that the number of observations n is inde¬ 
pendent from the chirp mass values we observe, we can 
rewrite this as 


P(U|{Al},n) 


P{{M}\n,V,)P{n\Vi)P{Vi) 

P{{M},n) 


( 20 ) 


We normalise by assuming that the discrete model vari¬ 
ations we consider cover all possible states of the uni¬ 
verse, which is an idealisation that we shall discuss in 
more detail later. However, this assumption allows us to 
define the normalisation factor by requiring the sum of 
the probabilities for each model to be unity, which leads 
to 


P(U|{M},n) 


P({M}|n,U)P(n|U)P(U) 

EfcU({M}|n,Ufe)P(n|Ufc)m)' ^ ^ 


We assume a uniform prior on the models, 

m) = ( 22 ) 

where Af is the number of models we are considering. 
The prior then cancels out and we are left with 


P{V,\{M},n) 


P{{M}\n,Vi)P{n\Vi) 

j:,Pi{M}\n,Vk)Pin\Vk)- 


(23) 


The likelihood of making n observations in a set time, 
given a model predicting mean number of observations, 
/Ti, is given by the Poisson distribution: 


P(n|U) = P{n\^ii) = (24) 


The likelihood terms of the form P({Al}|n, Vi) are cal¬ 
culated by binning the chirp mass distributions for each 
model into a histogram. We then calculate the likelihood 
of the observed samples being drawn from their bins us¬ 
ing the multinomial distribution 


O X}; 

P{{M}\n,V) = n\J[^-^ (25) 

1 1 
k—1 
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Figure 4. Chirp-mass distributions for each model in|Dominik et al.|(|2012|l using either their optimistic (top panel) or pessimistic (bottom 
pa nel) submodel and a 5 0-50 split of solar and sub-solar metallicities. me solid (blue online) line shows the intrinsic distribution as given 
by Dominik et al. ]( |2012| l. The dotted (green) line shows the same distribution when accounting for the observational biases introduced 
in becton [b.ij| as predicted for the early configuration of aLIGO. Finally, the dashed (red) line with largest chirp-mass support shows the 
expected observed distribution that additionally folds in the errors in measuring the chirp masses through GW observations. 
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where n is the number of samples in the observations, b 
is the number of bins, pik is the probability in model i 
of drawing a sample from bin k and Xk is the number of 
observations that fall into bin k, with 

'^Xk=n and = 1. (26) 

k k 


We calculate pik for each model and bin as the fraction 
of the total number of samples in the model which fall 
into that bi n. The bin size we employ is motived by 
Scott’s rule ( Scott|[T979 l, 


A = 


3.5cr 


(27) 


where A denotes the bin width, a is the standard devia¬ 
tion of the model, and is the total number of samples 
in model. To be able to consistently c omp are our simu¬ 
lated data with all models, we apply (27) to all models 
and then use the median bin width for the actual analy¬ 
sis. However, we find that changing this bin width by a 
factor of a few does not impact our results noticeably. 


5. OBSERVING SCENARIOS 


Table 2 

Predicted merger and detection rates. 


Vi 

TZ . 

B 

a 

A 

B A 

B 

(Ol) 

A 

fj ,’^ 1 

B 

(02) 

A 

0 

7.8 

40.8 

26.0 

24.9 

4.0 

25.2 

64 

402 

1 

4.6 

6.8 

27.3 

26.2 

2.3 

3.9 

37 

63 

2 

8.3 

36.0 

26.6 

24.9 

4.2 

25.9 

67 

413 

3 

4.0 

47.6 

25.0 

24.4 

1.9 

28.7 

30 

458 

4 

0.1 

3.1 

25.0 

24.7 

0.1 

1.9 

1 

30 

5 

7.8 

40.9 

26.0 

24.9 

4.0 

25.3 

64 

404 

6 

7.9 

41.3 

25.6 

24.2 

3.9 

25.1 

63 

401 

7 

8.6 

47.1 

25.3 

23.8 

4.0 

26.3 

65 

420 

8 

0.4 

2.1 

21.3 

10.0 

0.0 

0.6 

1 

9 

9 

11.8 

54.6 

23.2 

20.7 

3.4 

20.2 

54 

324 

10 

5.8 

31.3 

26.8 

26.2 

4.3 

26.0 

68 

415 

11 

10.4 

54.5 

29.8 

28.6 

8.5 

46.5 

136 

742 


Note. — The bina ry populations models, Vb predicted by 
|Dominik et al.|(|2012|l are summarised in Table (l] and the sub- 
models hi and A refer to pessimistic and optimistic assump¬ 
tions about the comm on envelope evolution of Hertzsprung 
gap donors (Sec. 

Local merger rate density in MWEG“^Myr“^. 

Average predicted observed chirp mass in Mq (see Sec.^ 

Mean number of detections predicted by each model for The 
early aLIGO observing runs Of and 02 (see text for details). 


The method we have developed transforms predicted 
binary distributions and merger rates into observable dis¬ 
tributions and detection rates which in turn can be con¬ 
fronted with a set of observations in order to assign pos¬ 
terior probabilities to each model. As such, the method is 
generally applicable to any set of theoretical predictions 
and detector configuration. 

In the following, however, we present results for specific 
choices of binary population models , det ector sensitivity 
and observing time. As detailed Sec. |2.2| and summarised 
in Table [T] we consid er 12 binary population models by 
Dominik et al. (2012), each with both the “pessimistic” 
(submodel H) and “optimistic” (submodel A) assump¬ 
tion about the common envelope evolution. This leads 
to 24 distinct predictions of the BBH chirp mass distri¬ 
bution (see Figure]^, where each comes with a distinct 
average merger rate density that we take from the arith¬ 
metic m ean of the solar and s ubsolar metallicity predic¬ 
tions by Dominik et al.| (2012) (Tables 2 and 3 therein). 
The loca!l merger rate densities for each model are given 
in Table Interestingly, due to the mass-dependent ob- 
servatioiial bias, models with higher merger rate density 
do not necessarily exhibit a higher detection rate, see for 
instance models 9 and 10 in Table [H 


Recent calculations by Dominik et al. (2013) that in¬ 
clude the cosmological evolution of merger rates give 
lower rate densities than the ones we infer from earlier 
work of the same authors. Consequently, the detection 
rates we find are up to a factor of 2 larger than those re¬ 


cently predicted by Dominik et al. 


(2014) (this is based 


on a direct comparison of our method with their other¬ 
wise equivalent approach using the same detector config¬ 
uration). However, this neither affects the general proof 
of principle carried out here, nor do the conclusions we 
shall draw in the following section change qualitatively 
by varying the detection rate at this level. 

We also have to specify in the sensitivity (i.e., noise 
spectral density) of our assumed GW detector and the 
observing time. Closely following Aasi et al. (|2013b ), we 
consider the first two aLIGO science runs dubbed U1 and 


02, respectively. The first science run for aLIGO (01) 
is planned to begin in autumn 2015. The duration of 01 
will be approximately 3 months for the two aLIGO de¬ 
tectors. We assume each detector has a duty cycle of 0.8 
so that the total period of coincident observation during 
01 will be about 0.16 years. The noise power s pectral 


density we u se is the “early aLIGO” conhguration (Shoe¬ 


maker 


2010 ). 


We further consider a second science run, 02. During 
02, the detectors are planned to observe for approxi¬ 
mately 6 months with a comparable duty cycle to 01. 
It is expected that, after further improvements of the 
instruments following 01, the aLIGO detectors during 
02 will be approximately a factor of 2 more sensitive 
than the nominal early aLIGO noise curve we use for 01. 
While the evolution of the noise power spectral density 
is in general a function of the frequency, we find that, 
in practic e, the difference be tween the predicted noise 
curves in Aasi et al. (2013b) results in improved hori¬ 


zon distaiices and error estimates that are well approx¬ 
imated by simply scaling the results we obtain for the 
early aLIGO configuration. Hence, we simulate 02 by 
multiplying the 01 horizon distance by 2. The Fisher- 
matrix errors change only due increased SNR at fixed 
distance. This increase in sensitivity leads to a factor of 
8 increase in volume meaning that, in total, 02 surveys 
16 times the time-volume of 01. We show in the follow¬ 
ing section that this is when we will begin to distinguish 
between astrophysical models. 

6. RESULTS: DISTINGUISHING BBH FORMATION 
MODELS 

6.1. First aLIGO observing run (01) 

We simulate the observed BBH system s, assuming the 
universe matches one of the models from IDominik et al.l 


(2012), and calculate the posterior probability for each 
model. We repeat the experiment 10000 times before 
turning to the next model to simulate the universe. Fig¬ 
ure gives the median posterior probability for each 
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Model 



Figu re 5. The median pos terior probability for each model in the 
set of|Dominik et aL|(|2012|| models after an 01 like observing pe¬ 
riod oi U.lt) years, calculated from 10000 repeats. The model which 
observations were drawn from is shown on the axis labelled Uni¬ 
verse. The models which these observations were then compared 
to is labelled Model, so that the probabilities in each row sum to 
one. Models 0-11 are described in Tab le 1. The two submodels, A 
and B, are described in Section|2.2.2| 


model. 

In cases where one or few models have a high prob¬ 
ability, these would be distinguishable from the other 
models. However, all models with a high probability 
would be consistent with the observations. We reiterate 
that here we restrict attention to the models in IDominikI 


et al. (2012). Of course these do not cover the full space 
of binary merger predictions. If we were to include a 
broader range of models, it is likely that the conclusions 
we are able to draw would be weaker as various models 
would lead to comparable rates and mass distributions. 
Nonetheless, some of the conclusions we reach, such as 
excluding a number of models if there are no observations 
in 01, are robust. 

We first observe that, for the most part, we would 
be able to distinguish between submodels A and B that 
corr espond to different common envelope scenarios (see 
Sec. 2.2.21. This is unsurprising as the predicted rates 
for the niajority of models are significantly higher for 
submodel A (cf. Table J^. Models which predict low de¬ 
tection rates for modeFA. remain degenerate with those 
in model B. The mass distribution from such a small 
sample does not provide enough additional information 
to break these degeneracies in the rates. For example, 
model 1 A uses a very high, fixed envelope binding en¬ 
ergy, meaning that most binaries entering a common en¬ 
velope event fail to throw off the common envelope and 
merge, causing them to never form BBH systems (for 


a mor e detailed discussion of this, see Dominik et al. 
(2012)). On the other hand, submodel B does not allow 
a binary to survive a common envelope event if the donor 
is on the Hertzsprung Gap, and so again, many binaries 
merge and never form BBHs. This generically lowers the 
merger rates and thus detection rates for submodel B 
models, leading to the degeneracy visible in the upper 
right quadrant of Fig. 


Another interesting example involves models 4 and 8 
that, in the pessimistic submodel B, are consistent with 
no observations at all during 01. Hence, they cannot 
be distinguished from each other, or indeed model 8 A, 
although they are favoured over all other models if indeed 
no detection are made. 

Within the two submodels, it is difficult to identify 
the correct model. Indeed, there are numerous varia¬ 
tions which would be indistinguishable from the standard 
model. The only model which can be clearly identified is 
model 11, a model which reduces the strength of stellar 
winds by a factor of 2 over the standard model. We now 
discuss why we are able to distinguish this model from 
the others in such a short observational period. 


6.2. Stellar winds 

In massive 0-type stars, stellar winds of high tem¬ 
perature charged gas are driven by radiation pressure. 
In Wolf-Rave t stars mass loss rates can be as high as 
lO^'^MQyr”^ (Nugis & Tamers 2002). This can cause 
stars to lose a large amount of mass prior to the super¬ 
nova. Theoretical uncertainties in modelling these mass 
loss rates therefore translate into uncertainties in the 
pre-su pernova masses for massive stars. |Dominik et al.| 
(2012) examine the effects of reducing the strength of 
stellar winds by a factor of 2 on the distribution of BBHs 
in their Variation 11. Firstly, reducing stellar winds re¬ 
sults in stars having a higher mass prior to supernova 
than they would otherwise have. This in turn leads to 
more mass falling back onto the compact object during 
formation, which reduces the magnitude of natal kicks 
given to black holes. This results in more systems sur¬ 
viving the supernova (rather than being disrupted) and 
increases the merger rates. More massive pre-supernova 
stars also form more massive remnants, resulting in the 
most massive BBH having a chirp mass of ~ 64M0 with 
reduced stellar winds compared to ~ STMq using the 
standard prescription. Finally, reducing the strength of 
stellar winds allows stars with a lower zero age main se¬ 
quence mass to form black holes due to more mass being 
retained. This can boost the BBH merger rate compared 
to the standard model. 

All of these effects combined mean that Variation 
11 predicts BBHs with characteristically higher chirp 
masses, as well as predicting a much higher merger rate 
than all other models (even for the pessimistic submodel 
B in 01, Variation 11 predicts 0(10) observations). We 
therefore expect that we would be able to correctly dis¬ 
tinguish a universe following Variation 11 from all other 
models with relatively few observations. In Figure we 
show the median posterior probability for each model 
as a function of the observation time, based on 10000 
redraws of the observations. We find that when draw¬ 
ing observations from a universe following Variation 11 
we overwhelmingly favour it within the duration of 01, 
with 0(10) observations. 


6.3. Second aLIGO observing run (02) 

We now turn our attention to the second observing run, 
02, and investigate which models can be distinguished 
using the much larger time-volume surveyed by 02. In 
Figure]^ we again show a matrix plot showing the (me¬ 
dian) posterior probability for each model after a period 
corresponding to the 02 run. 
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Figure 6. The median posterior probability for each of the models 
in the set as a function of observation time for a period of time cor¬ 
responding to the aLIGO 01 run (0.16 years). GW observations are 
drawn from a universe following Variation 11, submodel B which 
reduces the strength of stellar winds by a factor of 2 compared to 
the standard model. The blue (solid) line shows the median pos¬ 
terior probability for Variation 11 taken from 10000 repeats, and 
the shaded error bar shows the 68% confidence interval. Variations 
0,2,5,6,7 & 10 are plotted in green (dot-dash), while variations 1,3 
&: 9 are plotted in black (dotted). Variations 4 8 predicting ~ 0 

observations in 01 are plotted in red (dashed). 
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Figure 7 . The median poste rior probability for each model in 
the set of|Dominik et al.| l|2012|) models after an 02 like observing 
period of U.a2 years witn a detector more sensitive than the early 
aLIGO noise curve by a factor of 2. The median is calculated based 
on 10000 redraws of the observations. The model which observa¬ 
tions were drawn from is shown on the axis labelled Universe. The 
model which these observations were then compared to is labelled 
Model. Models 0-11 are described in Table[^ The two submodels, 
A and B, are described in Section|2.2.2| 


Figure has a more diagonal form than Figure [5 
meaning that in many cases the correct model is favourec 
and others are disfavoured within the 02 period. In par¬ 
ticular, the optimistic and pessimistic submodels A and 
B become almost enti rely distinct from each other. This 
is because most of the Dominik et al. (2012) models pre¬ 


dict 0(100) (0(10)) observations during the 02 period 
for the optimistic (pessimistic) submodels respectively 
(as shown in Table 2). Furthermore, the majority of 
variations in submodel A can be unambiguously identi¬ 
fied; the exception being that the standard model which 
remains degenerate with mode ls 5, 6 and 7, as we dis¬ 
cuss in detail in Section [6.3. 1[ For the pessimistic sub¬ 
model B, the standard model remains indistinguishable 
from a number of other variations. However, there are a 
few models which can be clearly distinguished, including 
models 4 and 8 (that predict significantly lower rates), 
and 9, 10 and 11. All of these models predict tens of 
observations and consequently, we are able to use infor¬ 
mation from both the chirp mass distribution and the 
detection rate to help distinguish models. Model 10 in¬ 
volves the variation of the su pernova engine, which we 
elaborate on in Section 16.3.21 

6.3.1. Black hole kicks and maximum neutron star mass 

Not all models are distinguishable, even with the 
0(100) observations predicted by the optimistic sub¬ 
model A for 02. For example, in Figurej^we see that the 
standard model is degenerate with Variations 5, 6 and 7. 
We now explain why this is so. 

As already mentioned, it is unclear what the correct 
distribution of natal kicks given to black holes upon for¬ 
mation is. In orde r to investigate the possibilities, Do¬ 
minik et al. (2012) vary two parameters relating to the 
kicks imparted onto newly formed black holes; the char¬ 
acteristic velocity a and the fraction of mass fh which 
falls back onto the newly born black hole. 

In their standard model, black holes receive a kick Vk 
whose magnitude Umax is drawn from a Maxwellian distri¬ 
bution with (T = 265km s~^, and reduced by the fraction 
of mass falling back onto the black hole ft as 

Vk = Umax(l - fb), (28) 

where ft is calcu lated using the prescription given in 


Fryer et al. 

(2012 

). 

In order to tes 

the effects of smaller natal kicks, in 

Variation 7 

Dominik et al. 

(2012 

) reduce the magnitude 


a factor of 2. They use a Maxwellian distribution with 
a = 132.5km s~^. For BBHs, this has very little effect 
on the chirp mass distribution, and so one cannot expect 
to be able to distinguish this model from one using full 
kicks. 

The same holds true when the maximum neutron-star 
mass is increased (decreased) from its fiducial value in 
the standard model of 2.5M0. This has very little im¬ 
pact on the BBH chirp mass distribution and so there 
is effectively a degeneracy between these models. This 
could be resolved by also including BNS observations in 
the comparison. We do not do this here as we concen- 
trate on the BBH pr edictions, due to the prediction by 


Dominik et al. (2012) that these will dominate the early 
aLlGU detections. 

6.3.2. Supernova engine 



mation, and thus the black hole masses. In particular. 
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they use the rapid supernova engine. When employed 
in a compact binary population code such as StarTrack, 
the r apid supernova engine reproduce s the observed mass 
gap (|Ozel et al.|2010t|Farr et al.|2011|) in compact objects 
between the highest mass neutron stars and the lowest 
mass black holes (for a discussion of using GW observa- 
tions to infer the pres e nce or absence of a mas s gap, see 


Hannam et al. (2013); Littenberg et al. (2015); Mandel 
et al.| [2015|). 


In model 10 Dominik et al. (2012) vary th is prescrip¬ 
tion to use the delayed supernova engine from [Fryer et al.| 
(2012), which produces a continuous distribution of black 
hole masses (and thus BBH chirp masses). We there¬ 
fore expect that the difference between these two models 
might be visible in the chirp mass distributions. We see 
however from Table that these two models predict sim¬ 
ilar merger rates for BBH, and so we do not expect to 
be able to distinguish them based on the detection rate. 
Nonetheless, we see from Figure that this model can 
be distinguished from the others oy the end of 02 and 
even, to a lesser degree, at the end of 01 (Figure 
To illustrate the importance of both the mass aistri- 
bution and predicted rates, in Figure!^ we show the re¬ 
sults that would be obtained using only one of these to 
separate the models. By comparing these results with 
Figure it becomes clear that both the mass and rate 
measurements contribute significantly to our ability to 
distinguish between models. As expected, the delayed 
supernova engine (model 10) is distinguished from ob¬ 
served masses, but the rates are quite degenerate with 
other models. In contrast, models 4 and 8, are distin¬ 
guished primarily by rate measurements, and not masses. 
As we have mentioned previously, the unknown spin dis¬ 
tribution of black holes in binary systems can change the 
rate by a factor of two or three. Similarly, both the mass 
and rate distributions are subject to uncertainties due 
to additional physical effects which are not yet incorpo¬ 
rated. Consequently, one might choose to incorporate an 
uncertainty in the rates or mass distributions. The re¬ 
sults in Figurel^illustrate the extreme scenario where one 
assumes knowledge of only the rate or mass distribution. 
Adding an uncertainty to the mass or rate distributions 
will lead to a result between those shown in Figure and 


7. SUMMARY AND FUTURE WORK 

In this paper we have outlined a method for comparing 
GW observations of BBH mergers to binary population 
synthesis predictions using a Bayesian model comparison 
framewor k. Starting from chirp mass distribution pre¬ 
dicted by Dominik et al. (2012), we produce predicted 
observed chirp mass distributions accounting for known 
observational effects. We incorporate 


(a) The redshifting of observed binary masses due to the 
cosmological distances out to which they will be ob¬ 
served. 


(b) The observational bias of GW detectors to detect 
more massive systems, since they can be seen to 
greater distances and thus in much larger volumes. 

(c) Fisher matrix estimates of measurement uncertain¬ 
ties in the recovery of the chirp mass of BBH. 


We sho w that given the merger rates predicted by the 
models of Dominik et al. (2012), we will begin to be 
able to distinguish between population synthesis models 
in the first two aLlGO science runs. Ruling out models 
in turn can help to constrain the value of unknown pa¬ 
rameters, which relate to poorly understood astrophysics 
relating to binary evolution. 

Of course, the set of models considered here by no 
means encompasses the full set of stellar evolution models 
available in the literature. We restricted attention to this 
subset of models as the data was publicly available in an 
easy to use form. It would be straightforward to include 
additional models into this analysis. Ideally, we would 
make use of a dense set of models, where numerous astro- 
physical parameters are jointly varied. This would allow 
us to interp olate between models, and extract best-fit 
param eters (O’Shaughnessy 2013 O’Shaughnessy et al.| 


2008). Furthermore, we have restricted attention to the 
two best-measured quantities: the rate and chirp mass 
distribution of binaries, and only used point estimates 
of the masses. The inclusion of full parameter distribu¬ 
tions can only enhance our ability to distinguish between 
models. 

The method we have introduced allows us to distin¬ 
guish between a given set of stellar evolution models. It 
will identify the model, or models, that best agree with 
the observed rate and mass distribution. It will not, how¬ 
ever, indicate whether the best model is actually a good 
fit to the observations — only that it is better than the 
others. This could be remedied by introducing a simple, 
generic model. For example, the intrinsic mass distribu¬ 
tions shown in Figure are reasonably well desribed by 
a decaying power law with an upper and lower mass cut¬ 
off. One could then imagine extending the set of mod¬ 
els to include this phenomenological mass distribution 
parametrized by three variables with an additional vari¬ 
able rate. To calculate the posterior for this distribution, 
we would then have to marginalize over four parameters. 
Thus, even if the generic model was a reasonable fit to the 
data, it would be penalized by the large initial parameter 
space. It is likely that the generic model would be pre¬ 
ferred after a small number of observations. With a large 
number of observations, the rate and mass distributions 
would be reasonably well measured. Any specific model 
which matched the observations well would then be pre¬ 
ferred to the generic model due to its broader support on 
the parameter space. It would be reasonably straightfor¬ 
ward to extend our method to include a generic model, 
and this is something we plan to incorporate in the fu¬ 
ture. 

In this study we concentrated on the information that 
could be gained from GW observations of BBH mergers. 
aLIGO and AdV are also expected to observe the inspi¬ 
ral, merger and ringdown of compact binaries including 
neutron stars (BNS and NSBH systems). One should in¬ 
clude all GW observations of compact binaries in order 
to extract the maximum amount of information from the 
observations. In fact, as discussed above, we are unable 
to distinguish models which vary the maximum allowed 
neutron star mass since we ignore these events here. In 
this study we ignored these events since the predicted de¬ 
tection rates for BBH mergers dominated those of other 
compact binary mergers. The BBH mass distribution 
also spans a large range of masses, with structure encod- 
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Figure 8. Probabilities for the scenario of Fig.j^ separated into contributions from the rates (left) and the mass distribution (right). 
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ing information about binary evolution. Ignoring other 
families of compact binaries also allowed us to avoid am¬ 
biguities in discerning the family of the source (BNS, 
NSBH or BBH) due to degeneracies whic h exist in mea¬ 


surin g the mass ratio for these systems (Hannam et al. 
20131), although this can be dealt with in the future (|h'arr 


eFaf||20T3 ). 


All these considerations have to be carefully taken into 
account in future studies. However, our results indicate 
that the upcoming generation of advanced GW detectors 
will soon start putting non-trivial bounds on current and 
future binary evolution models, and analyses like the one 
presented here will provide an important basis to link 
theoretical models with GW observations. 
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